/* Copyright (C) 2009, Lukasz Kaczmarczyk (likask AT civil.gla.ac.uk)
 * --------------------------------------------------------------
 * FIXME: DESCRIPTION
 */

/* This file is part of yaffems.
 * yaffems is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * yaffems is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with yaffems. If not, see <http://www.gnu.org/licenses/>. */

#include <driver.hpp>
#include<math.h>
#include<complex>

PetscErrorCode ierr;
ErrorCode rval;
//int imesh_err;
PETSCTRUTH flg;

#define MYPCOMM_INDEX 0
void driver_GetFromOption(snes_ctx *ctx) {
  PETSCTRUTH flg_common = PETSC_TRUE;
  PetscOptionsGetString(PETSC_NULL,"-my_file_name",ctx->mesh_file_name,255,&flg);
  if(flg != PETSC_TRUE) {
    PetscPrintf(PETSC_COMM_WORLD,"*** ERROR -my_file_name (MESH FILE NEEDED)\n");
    flg_common = PETSC_FALSE; }
  //MATERIAL PARAMETERS
  double E;
  PetscOptionsGetReal(PETSC_NULL,"-my_E",&E,&flg);
  if( flg != PETSC_TRUE ) {
    PetscPrintf(PETSC_COMM_WORLD,"*** ERROR -my_E (YOUNG MODULUS NEEDED)\n");
    flg_common = PETSC_FALSE; 
  }
  double nu;
  PetscOptionsGetReal(PETSC_NULL,"-my_nu",&nu,&flg);
  if( flg != PETSC_TRUE ) {
    PetscPrintf(PETSC_COMM_WORLD,"*** ERROR -my_nu (POISON RATIO NEEDED)\n");
    flg_common = PETSC_FALSE; 
  }
  ctx->E = E;
  ctx->nu = nu;
  ctx->lambda = LAMBDA(E,nu);
  ctx->mu = MU(E,nu);
  PetscOptionsGetReal(PETSC_NULL,"-my_eps",&ctx->eps,&flg);
  if( flg == PETSC_FALSE ) ctx->eps = 1e-12;
  PetscOptionsGetReal(PETSC_NULL,"-my_gc",&ctx->gc,&flg);
  if( flg == PETSC_FALSE ) ctx->gc;
  PetscOptionsGetReal(PETSC_NULL,"-my_gc_threshold",&ctx->gc_threshold,&flg);
  if( flg == PETSC_FALSE ) ctx->gc_threshold = 0.5;
  //PetscOptionsGetReal(PETSC_NULL,"-my_crack_barrier0",&ctx->penalty_barrer0,&flg);
  //if( flg == PETSC_FALSE ) ctx->penalty_barrer0 = 0;
  //PetscOptionsGetReal(PETSC_NULL,"-my_crack_barrier1",&ctx->penalty_barrer1,&flg);
  if( flg == PETSC_FALSE ) ctx->penalty_barrer1 = 0;
  //MESH ADAPATION PARAMETERS
  PetscOptionsGetReal(PETSC_NULL,"-my_gamma",&ctx->gamma,&flg);
  if( flg == PETSC_FALSE ) ctx->gamma = 0.1;
  PetscOptionsGetReal(PETSC_NULL,"-my_alpha",&ctx->alpha,&flg);
  if( flg == PETSC_FALSE ) ctx->alpha = 0;
  PetscOptionsGetReal(PETSC_NULL,"-my_alpha_min",&ctx->alpha_min,&flg);
  if( flg == PETSC_FALSE ) ctx->alpha_min = ctx->alpha;
  PetscOptionsGetReal(PETSC_NULL,"-my_beta",&ctx->beta,&flg);
  if( flg == PETSC_FALSE ) ctx->beta = 0.5;
  PetscOptionsGetReal(PETSC_NULL,"-my_kappa",&ctx->kappa,&flg);
  if( flg == PETSC_FALSE ) ctx->kappa = 1;
  PetscOptionsGetReal(PETSC_NULL,"-my_cred",&ctx->c_red,&flg);
  if( flg == PETSC_FALSE ) ctx->c_red = 0.6;
  PetscOptionsGetReal(PETSC_NULL,"-my_cinc",&ctx->c_inc,&flg);
  if( flg == PETSC_FALSE ) ctx->c_inc = 1./ctx->c_red;
  assert( flg_common == PETSC_TRUE );
  ctx->moab = NULL;
  ctx->meshset = 0;
  ctx->is_hierarchical = 0;
  ctx->joker_bc = FREE_TO_MOVE;
  ctx->arc_diagonal = 0;
  ctx->arc_length_value = 0;
  bzero(ctx->arc_length_direction,6*sizeof(double));
}
void driver_partition(Interface& moab) {
  make_clear_all_data();
  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
  if(pcomm == NULL) pcomm =  new ParallelComm(&moab,PETSC_COMM_WORLD);
  Tag ptag;
  int def_val = -1;
  //rval = moab.tag_get_handle(PARALLEL_PARTITION_TAG_NAME,ptag); CHKERR(rval);
  moab.tag_create(PARALLEL_PARTITION_TAG_NAME,sizeof(int),MB_TAG_DENSE,MB_TYPE_INTEGER,ptag,&def_val,true);
  Tag th_poly_order;
  int default_int = 1;
  moab.tag_create(POLY_ORDER,sizeof(int),MB_TAG_SPARSE,MB_TYPE_INTEGER,th_poly_order,&default_int,true); 
  Tag th_moab_gid;
  moab.tag_create(PARALLEL_GID_TAG_NAME,sizeof(int),MB_TAG_DENSE,MB_TYPE_INTEGER,th_moab_gid,&def_val,true);
  make_aux_moab_ents_numeration_by_type(moab,0,MBTET,STANDARD_SOLID_FE);
  Range elems;
  //moab.get_entities_by_type(0,MBTET,elems); 
  make_aux_get_entities_by_ents_type(moab,0,MBTET,STANDARD_SOLID_FE,elems); 
  //moab.get_entities_by_dimension(0,3,elems); 
  int local_nb = (int)ceil((double)elems.size()/pcomm->size());
  vector<int> i,j;
  for(Range::iterator eit = elems.begin();eit!=elems.end();eit++) {
    int elem_gid;
    moab.tag_get_data(th_moab_gid,&*eit,1,&elem_gid); 
    elem_gid--;
    if(elem_gid<((int)pcomm->rank()*local_nb)) continue;
    if(elem_gid>=((int)(pcomm->rank()+1)*local_nb)) continue;
    i.push_back(j.size());
    Range nodes; 
    moab.get_connectivity(&*eit,1,nodes,true); 
    Range adj_elems;
    make_aux_get_adjacencies(moab,0,MBTET,STANDARD_SOLID_FE,nodes,3,true,adj_elems,Interface::UNION);
    adj_elems.erase(*eit);
    Range::iterator eit2 = adj_elems.begin();
    for(;eit2!=adj_elems.end();eit2++) {
      int elem2_gid;
      moab.tag_get_data(th_moab_gid,&*eit2,1,&elem2_gid); 
      elem2_gid--;
      assert(elem2_gid!=elem_gid);
      j.push_back(elem2_gid);
    }
  }
  i.push_back(j.size());
  //fprintf(stderr,"elems.size() %d i.size() %d j.size() %d\n",(int)elems.size(),(int)i.size(),(int)j.size());
  //PetscPrintf(PETSC_COMM_WORLD,"elems.size() %d i.size() %d j.size() %d\n",(int)elems.size(),(int)i.size(),(int)j.size());
  //fprintf(stderr,"%d %d %d\n",i[3227-1],i[3227],i[3227+1]);
  int *ii,*jj;
  PetscMalloc(sizeof(int)*i.size(),&ii);
  bcopy(&i[0],ii,sizeof(int)*i.size());
  PetscMalloc(sizeof(int)*j.size(),&jj);
  bcopy(&j[0],jj,sizeof(int)*j.size()); 
  Mat Adj;
  MatCreateMPIAdj(pcomm->comm(),i.size()-1,elems.size(),ii,jj,PETSC_NULL,&Adj);
  MatPartitioning part;
  MatPartitioningCreate(pcomm->comm(),&part);
  MatPartitioningSetAdjacency(part,Adj);
  int *elems_poly_order;
  PetscMalloc(sizeof(int)*elems.size(),&elems_poly_order);
  rval = moab.tag_get_data(th_poly_order,elems,elems_poly_order); CHKERR(rval);
  for(int ee = 0;ee<(int)elems.size();ee++) {
    elems_poly_order[ee] = 4 + 6*NBEDGE_H1(elems_poly_order[ee]) + 4*NBFACE_H1(elems_poly_order[ee]) + NBVOLUME_H1(elems_poly_order[ee]);
  }
  MatPartitioningSetVertexWeights(part,elems_poly_order);
  MatPartitioningSetFromOptions(part);
  IS is;
  MatPartitioningApply(part,&is);
  //PetscPrintf(PETSC_COMM_WORLD,"New processor assignment for each element\n");
  //ISView(is,PETSC_VIEWER_STDOUT_WORLD);
  MATPARTITIONINGDESTROY(part);
  MATDESTROY(Adj);
  int count[pcomm->size()];
  ISPartitioningCount(is,pcomm->size(),count);
  //IS isg;
  //ISPartitioningToNumbering(is,&isg);
  IS is_gather;
  ISAllGather(is,&is_gather);
  const int *part_number;
  ISGetIndices(is_gather,&part_number); 
  rval = moab.tag_set_data(ptag,elems,part_number); CHKERR(rval);
  ISRestoreIndices(is_gather,&part_number); 
  ISDESTROY(is);
  //ISDESTROY(isg);
  ISDESTROY(is_gather);
}
#define DRIVER_SETPSTAT(CONDITION) \
    part = rr; \
    Range notown_elems; \
    Range::iterator elem_it = all_elems.begin(); \
    for(;elem_it!=all_elems.end();elem_it++) { \
      int elem_part; \
      moab.tag_get_data(ptag,&*elem_it,1,&elem_part); \
      if(elem_part == part) notown_elems.insert(*elem_it); \
    } \
    elem_it = notown_elems.begin(); \
    for(;elem_it!=notown_elems.end();elem_it++) { \
      unsigned int pstat = PSTATUS_NOT_OWNED; \
      moab.tag_set_data(pstatusTag,&*elem_it,1,&pstat); \
    } \
    Range notown_nodes; \
    moab.get_connectivity(notown_elems,notown_nodes,true); \
    if(CONDITION) notown_nodes = subtract(notown_nodes,nodes); \
    Range::iterator nit = notown_nodes.begin(); \
    for(;nit!=notown_nodes.end();nit++) { \
      unsigned int pstat = PSTATUS_NOT_OWNED; \
      moab.tag_set_data(pstatusTag,&*nit,1,&pstat); \
    } \
    Range notown_faces; \
    make_aux_get_adjacencies(moab,0,MBTRI,STANDARD_SOLID_FE_BRIDGE,notown_elems,2,true,notown_faces,Interface::UNION); \
    if(CONDITION) notown_faces = subtract(notown_faces,faces); \
    Range::iterator fit = notown_faces.begin(); \
    for(;fit!=notown_faces.end();fit++) { \
      unsigned int pstat = PSTATUS_NOT_OWNED; \
      moab.tag_set_data(pstatusTag,&*fit,1,&pstat); \
    } \
    Range notown_edges; \
    make_aux_get_adjacencies(moab,0,MBEDGE,STANDARD_SOLID_FE_BRIDGE,notown_elems,1,true,notown_edges,Interface::UNION); \
    if(CONDITION) notown_edges = subtract(notown_edges,edges); \
    Range::iterator eit = notown_edges.begin(); \
    for(;eit!=notown_edges.end();eit++) { \
      unsigned int pstat = PSTATUS_NOT_OWNED; \
      moab.tag_set_data(pstatusTag,&*eit,1,&pstat); \
    }
void driver_create_meshset(Interface& moab,EntityHandle meshset,snes_ctx *ctx) {
  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
  if(pcomm == NULL) pcomm =  new ParallelComm(&moab,PETSC_COMM_WORLD);
  moab.clear_meshset(&meshset,1);
  make_clear_all_data();
  make_aux_set_bridge_type_by_type(moab,0,MBTRI,2,STANDARD_SOLID_FE,STANDARD_SOLID_FE_BRIDGE);
  make_aux_set_bridge_type_by_type(moab,0,MBEDGE,1,STANDARD_SOLID_FE,STANDARD_SOLID_FE_BRIDGE);
  Range elems,faces,edges,nodes;
  moab.get_entities_by_dimension(0,3,elems,true); 
  if(0/*clening_mesh*/) {
    //cleaning mesh
    make_aux_get_range_subset_by_type(moab,(EntityType)-1,STANDARD_SOLID_FE|INTERFACE_SOLID_FE,elems);
    make_aux_get_adjacencies(moab,0,(EntityType)(-1),STANDARD_SOLID_FE_BRIDGE|INTERFACE_SOLID_FE_BRIDGE,elems,2,true,faces,Interface::UNION); 
    make_aux_get_adjacencies(moab,0,(EntityType)(-1),STANDARD_SOLID_FE_BRIDGE|INTERFACE_SOLID_FE_BRIDGE,elems,1,true,edges,Interface::UNION); 
    Range edges_all_subtract;
    moab.get_entities_by_dimension(0,1,edges_all_subtract,true); 
    Range faces_all_subtract;
    moab.get_entities_by_dimension(0,2,faces_all_subtract,true); 
    edges_all_subtract = subtract(edges_all_subtract,edges);
    if(edges_all_subtract.size() !=0 ) {
      PetscPrintf(pcomm->comm(),"**** WARRNING: hanging edges %u\n",edges_all_subtract.size());
      /*for(Range::iterator eit = edges_all_subtract.begin();eit!=edges_all_subtract.end();eit++) {
        const EntityHandle* conn; 
        int num_nodes; 
        moab.get_connectivity(*eit,conn,num_nodes,true); 
        assert(num_nodes==2);
        fprintf(stderr,"ee %d: %d %d\n",(int)*eit,(int)conn[0],(int)conn[1]);
      }*/
      rval = moab.delete_entities(edges_all_subtract);  CHKERR(rval);
      make_clear_all_data();
    }
    //assert(edges_all_subtract.size()==0);
    faces_all_subtract = subtract(faces_all_subtract,faces);
    if(faces_all_subtract.size() !=0 ) {
      PetscPrintf(pcomm->comm(),"**** WARRNING: hanging faces %u\n",faces_all_subtract.size());
      //for(Range::iterator fit = faces_all_subtract.begin();fit!=faces_all_subtract.end();fit++) {
       // fprintf(stderr,"ff %d\n",(int)*fit);
      //}
      rval = moab.delete_entities(faces_all_subtract);  CHKERR(rval);
      make_clear_all_data();
    }
    rval = moab.get_connectivity(elems,nodes,true); CHKERR(rval); 
    Range all_nodes_subtracted;
    moab.get_entities_by_type(0,MBVERTEX,all_nodes_subtracted); 
    all_nodes_subtracted = subtract(all_nodes_subtracted,nodes);
    if(all_nodes_subtracted.size() != 0) {
      PetscPrintf(pcomm->comm(),"**** WARRNING: hanging nodes  %u\n",all_nodes_subtracted.size());
      for(Range::iterator nit = all_nodes_subtracted.begin();nit!=all_nodes_subtracted.end();nit++) {
        make_print_bc_tags(stderr,moab,*nit);
      }
      rval = moab.delete_entities(all_nodes_subtracted);  CHKERR(rval);
      make_clear_all_data();
    }
    //assert(faces_all_subtract.size()==0);
    elems.clear();
    faces.clear();
    edges.clear();
    nodes.clear();
    //check mesh <- very time consiuming
    //Core *core = (Core *)(&moab);
    //rval = core->check_adjacencies(); CHKERR(rval);
    //cleaning mesh finish
  }
  Tag ptag;
  rval = moab.tag_get_handle(PARALLEL_PARTITION_TAG_NAME,ptag); CHKERR(rval);
  //set part tags to interface elements
  Range int_elems;
  moab.get_entities_by_dimension(0,3,int_elems,true); 
  make_aux_get_range_subset_by_type(moab,(EntityType)-1,INTERFACE_SOLID_FE,int_elems);
  for(Range::iterator eit = int_elems.begin();eit!=int_elems.end();eit++) {
    Range nodes; 
    moab.get_connectivity(&*eit,1,nodes,true); 
    Range adj_elems;
    make_aux_get_adjacencies(moab,0,MBTET,STANDARD_SOLID_FE,nodes,3,true,adj_elems,Interface::UNION);
    int part_number = 0;
    for(Range::iterator eit2 = adj_elems.begin();eit2!=adj_elems.end();eit2++) {
      int part_;
      rval = moab.tag_get_data(ptag,&*eit2,1,&part_); CHKERR(rval);
      part_number = part_number > part_ ? part_number : part_;
    }
    rval = moab.tag_set_data(ptag,&*eit,1,&part_number); CHKERR(rval);
  }
  //finish set part tags to interface elements
  Tag tag_handles[] = { ptag };
  int part = pcomm->rank();
  void* tag_data[] = { (void*)&part };
  moab.get_entities_by_type_and_tag(0,MBTET,tag_handles,tag_data,1,elems);
  moab.add_entities(meshset,elems);
  make_aux_get_adjacencies(moab,0,(EntityType)(-1),STANDARD_SOLID_FE_BRIDGE,elems,2,true,faces,Interface::UNION); 
  moab.add_entities(meshset,faces);
  make_aux_get_adjacencies(moab,0,(EntityType)(-1),STANDARD_SOLID_FE_BRIDGE,elems,1,true,edges,Interface::UNION); 
  moab.add_entities(meshset,edges);
  rval = moab.get_connectivity(elems,nodes,true); CHKERR(rval); 
  moab.add_entities(meshset,nodes);
  make_aux_moab_ents_numeration_by_type(moab,0,MBVERTEX,NO_DEFINED_ELEM);
  make_aux_moab_ents_numeration_by_type(moab,0,MBTET,STANDARD_SOLID_FE); 
  Tag pstatusTag;
  unsigned int tmp_pstatus = 0;
  moab.tag_create(MY_PARALLEL_STATUS_TAG_NAME,sizeof(unsigned int),MB_TAG_DENSE,MB_TYPE_OPAQUE,pstatusTag,&tmp_pstatus,true);
  //clean pstat tags
  for(Range::iterator eit = elems.begin();eit!=elems.end();eit++) {
    unsigned int pstat = 0;
    moab.tag_set_data(pstatusTag,&*eit,1,&pstat); 
  }
  for(Range::iterator nit = nodes.begin();nit!=nodes.end();nit++) {
    unsigned int pstat = 0;
    moab.tag_set_data(pstatusTag,&*nit,1,&pstat); 
  }
  Range::iterator eit = faces.begin();
  for(;eit!=faces.end();eit++) {
    unsigned int pstat = 0;
    moab.tag_set_data(pstatusTag,&*eit,1,&pstat); 
  }
  eit = edges.begin();
  for(;eit!=edges.end();eit++) {
    unsigned int pstat = 0;
    moab.tag_set_data(pstatusTag,&*eit,1,&pstat); 
  }
  //set pstat tags
  unsigned int rr = 0;
  Range all_elems;
  moab.get_entities_by_dimension(0,3,all_elems,true); 
  for(;rr<pcomm->rank();rr++) {
    //those are own
    DRIVER_SETPSTAT(0) 
  }
  for(rr++;rr<pcomm->size();rr++) {
    //those are not own
    DRIVER_SETPSTAT(1) 
  }
  ctx->moab = &moab;
  ctx->meshset = meshset;
  make_aux_calulate_volume_and_quality(moab,meshset,&ctx->Vol,&ctx->ave_Vol,&ctx->min_quali);
  PetscPrintf(pcomm->comm(),"Vol = %6.4e ave_Vol = %6.4e min_quali = %6.4e\n",ctx->Vol,ctx->ave_Vol,ctx->min_quali);
}
PetscErrorCode driver_create_matrices(Interface& moab,EntityHandle meshset,Mat *M,Mat *B,Vec *F,Vec *D,Vec *Lambda,Vec *FfromArray,Vec *DfromArray,snes_ctx *ctx) {
  make_aux_petsc_node_numeration(moab,meshset);
  make_aux_petsc_elem_numeration(moab,meshset);

  int nb_lambda_dofs;
  make_aux_nb_Lambda_dofs(moab,meshset,&nb_lambda_dofs);
  make_aux_petsc_vecANDmatrix(moab,meshset,B,F,PETSC_NULL,FfromArray,nb_lambda_dofs);
  make_aux_petsc_vecANDmatrix(moab,meshset,PETSC_NULL,D,Lambda,DfromArray,nb_lambda_dofs);
  VecZeroEntries(*FfromArray);
  VecZeroEntries(*DfromArray);
  ctx->B = *B;
  ctx->Bprec = *B;
  snes_ksp_mat_shell(M,*FfromArray,ctx);
  ctx->is_hierarchical = 0;
  return 0;
}
PetscErrorCode driver_destroy_matrices(Mat M,Mat B,Vec F,Vec D,Vec Lambda,Vec FfromArray,Vec DfromArray,snes_ctx *ctx) {
  MATDESTROY(M);
  MATDESTROY(B);
  VECDESTROY(F);
  VECDESTROY(D);
  VECDESTROY(Lambda);
  VECDESTROY(FfromArray);
  VECDESTROY(DfromArray);
  return 0;
}
PetscErrorCode driver_create_matrices_hierarchical(Interface& moab,EntityHandle meshset,Mat *M,Mat *B,Vec *F,Vec *D,Vec *F_hierarchical,Vec *D_hierarchical,Vec *Lambda,Vec *FfromArray,Vec *DfromArray,snes_ctx *ctx) {
  make_aux_petsc_node_numeration(moab,meshset);
  make_aux_petsc_elem_numeration(moab,meshset);
  make_hierarchical_set_edge_poly_order_min(moab,meshset);
  make_hierarchical_set_face_poly_order_min(moab,meshset);
  make_hierarchical_moab_edge_numeration(moab,0);
  make_hierarchical_moab_face_numeration(moab,0);
  make_hierarchical_petsc_edge_numeration(moab,meshset);
  make_hierarchical_petsc_face_numeration(moab,meshset);
  make_hierarchical_petsc_edge_dofs_numeration(moab,meshset);
  make_hierarchical_petsc_face_dofs_numeration(moab,meshset);
  make_hierarchical_petsc_elem_dofs_numeration(moab,meshset);
  int nb_lambda_dofs;
  make_aux_nb_Lambda_dofs(moab,meshset,&nb_lambda_dofs);
  make_hierarchical_petsc_vecANDmatrix(moab,meshset,B,D_hierarchical,Lambda,DfromArray,nb_lambda_dofs);
  make_aux_petsc_createGhostVecArray_formLargeVec(moab,meshset,*DfromArray,D,nb_lambda_dofs);
  VecDuplicate(*DfromArray,FfromArray);
  make_aux_petsc_createGhostVecArray_formLargeVec(moab,meshset,*FfromArray,F,nb_lambda_dofs);
  make_hierarchical_petsc_createGhostVecArray_formLargeVec(moab,meshset,*FfromArray,F_hierarchical,nb_lambda_dofs);
  VecZeroEntries(*FfromArray);
  VecZeroEntries(*DfromArray);
  ctx->B = *B;
  ctx->Bprec = *B;
  snes_ksp_mat_shell(M,*FfromArray,ctx);
  ctx->is_hierarchical = 1;
  return 0;
}
PetscErrorCode driver_destroy_matrices_hierarchical(Mat M,Mat B,Vec F,Vec D,Vec F_hierarchical,Vec D_hierarchical,Vec Lambda,Vec FfromArray,Vec DfromArray,snes_ctx *ctx) {
  MATDESTROY(M);
  MATDESTROY(B);
  VECDESTROY(F);
  VECDESTROY(D);
  VECDESTROY(F_hierarchical);
  VECDESTROY(D_hierarchical);
  VECDESTROY(Lambda);
  VECDESTROY(FfromArray);
  VECDESTROY(DfromArray);
  return 0;
}
PetscErrorCode driver_setup_snes(Mat M,Mat B,Vec Dofs,Vec F,snes_ctx *ctx,SNES *snes) {
  ParallelComm* pcomm = ParallelComm::get_pcomm(ctx->moab,MYPCOMM_INDEX);
  ierr = SNESCreate(pcomm->comm(),snes); CHKERRQ(ierr);
  ierr = SNESSetApplicationContext(*snes,ctx); CHKERRQ(ierr);
  ierr = SNESSetFunction(*snes,F,rhs_snes,ctx); CHKERRQ(ierr);
  ierr = SNESSetJacobian(*snes,M,M,lhs_snes,ctx); CHKERRQ(ierr);
  ierr = SNESSetFromOptions(*snes); CHKERRQ(ierr);
  ierr = SNESSetConvergenceTest(*snes,snes_convernce_test,ctx,PETSC_NULL); CHKERRQ(ierr);
  ierr = SNESMonitorSet(*snes,snes_solution_monitor,ctx,PETSC_NULL); CHKERRQ(ierr);
  ctx->snes = *snes;
  ierr = SNESGetKSP(*snes,&(ctx->ksp)); CHKERRQ(ierr);
  ierr = KSPMonitorSet(ctx->ksp,ksp_soultion_monitor,ctx,PETSC_NULL); CHKERRQ(ierr);
  //create implicit preconditioner
  ierr = PCCreate(pcomm->comm(),&(ctx->pc)); CHKERRQ(ierr);
  ierr = PCCreate(pcomm->comm(),&(ctx->pc_uzawa)); CHKERRQ(ierr);
  #ifdef PETSC_3_3_
  //setup linesearch
  ierr = SNESLineSearchCreate(pcomm->comm(),&ctx->linesearch);  CHKERRQ(ierr);
  ierr = SNESLineSearchSetFromOptions(ctx->linesearch); CHKERRQ(ierr);
  ierr = SNESLineSearchSetSNES(ctx->linesearch,*snes); CHKERRQ(ierr);
  SNESLineSearch linesearch;
  ierr = SNESGetSNESLineSearch(*snes,&linesearch); CHKERRQ(ierr);
  ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHSHELL); CHKERRQ(ierr);
  ierr = SNESLineSearchShellSetUserFunc(linesearch,LineSearch,ctx); CHKERRQ(ierr);
  ierr = SNESSetSNESLineSearch(*snes,linesearch); CHKERRQ(ierr);
  #endif
  //setup my preconditiones
  PC pc;
  ierr = KSPGetPC(ctx->ksp,&pc); CHKERRQ(ierr);
  ierr = PCSetType(pc,PCSHELL); CHKERRQ(ierr);
  ierr = PCShellSetContext(pc,ctx); CHKERRQ(ierr);
  ierr = PCShellSetApply(pc,pc_uzawa_apply_arc_length); CHKERRQ(ierr);
  ierr = PCShellSetSetUp(pc,pc_uzawa_setup_arc_length); CHKERRQ(ierr);
  ierr = VecDuplicate(F,&(ctx->b)); CHKERRQ(ierr);
  ierr = VecDuplicate(ctx->b,&ctx->d_b); CHKERRQ(ierr);
  ierr = VecDuplicate(F,&(ctx->delta)); CHKERRQ(ierr);
  ierr = VecDuplicate(F,&(ctx->dF_dLambda)); CHKERRQ(ierr);
  ierr = VecDuplicate(F,&(ctx->dD_dF_dLambda)); CHKERRQ(ierr);
  VecZeroEntries(ctx->b);
  VecZeroEntries(ctx->d_b);
  VecZeroEntries(ctx->delta);
  VecZeroEntries(ctx->dF_dLambda);
  VecZeroEntries(ctx->dD_dF_dLambda);
  ierr = VecAssemblyBegin(ctx->delta); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(ctx->delta); CHKERRQ(ierr);
  return 0;
}
PetscErrorCode driver_setup_arc_length(Vec Dofs,snes_ctx *ctx) {
  ParallelComm* pcomm = ParallelComm::get_pcomm(ctx->moab,MYPCOMM_INDEX);
  PetscErrorCode ret_val = 0;
  VecZeroEntries(ctx->b);
  VecAssemblyBegin(ctx->b);
  VecAssemblyEnd (ctx->b);
  int nb_lambda_dofs;
  make_aux_nb_Lambda_dofs(*(ctx->moab),ctx->meshset,&nb_lambda_dofs);
  Vec Ghost_b,Ghost_b_hierarchical,GhostDofs;
  if(ctx->is_hierarchical) {
    make_hierarchical_petsc_createGhostVecArray_formLargeVec(*(ctx->moab),ctx->meshset,ctx->b,&Ghost_b_hierarchical,nb_lambda_dofs);
  }
  make_aux_petsc_createGhostVecArray_formLargeVec(*(ctx->moab),ctx->meshset,ctx->b,&Ghost_b,nb_lambda_dofs);
  make_aux_petsc_createGhostVecArray_formLargeVec(*(ctx->moab),ctx->meshset,Dofs,&GhostDofs,nb_lambda_dofs);
  if(get_arc_area_control()) {
    make_area_calculate_clear_data();
    make_area_calculate_values(ctx->eps,*(ctx->moab),ctx->meshset,GhostDofs,MATERIAL_M,true);
    if(ctx->is_hierarchical) {
      VecGhostUpdateBegin(Ghost_b_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
      VecGhostUpdateEnd(Ghost_b_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
      make_area_arc_length_b_vector(Ghost_b_hierarchical);
      VecGhostUpdateBegin(Ghost_b_hierarchical,ADD_VALUES,SCATTER_REVERSE);
      VecGhostUpdateEnd(Ghost_b_hierarchical,ADD_VALUES,SCATTER_REVERSE);
      VecGhostUpdateBegin(Ghost_b_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
      VecGhostUpdateEnd(Ghost_b_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
    } else {  
      VecGhostUpdateBegin(Ghost_b,INSERT_VALUES,SCATTER_FORWARD);
      VecGhostUpdateEnd(Ghost_b,INSERT_VALUES,SCATTER_FORWARD);
      make_area_arc_length_b_vector(Ghost_b);
      VecGhostUpdateBegin(Ghost_b,ADD_VALUES,SCATTER_REVERSE);
      VecGhostUpdateEnd(Ghost_b,ADD_VALUES,SCATTER_REVERSE);
      VecGhostUpdateBegin(Ghost_b,INSERT_VALUES,SCATTER_FORWARD);
      VecGhostUpdateEnd(Ghost_b,INSERT_VALUES,SCATTER_FORWARD);
    }
    double nrm2_dofs,nrm2_b;
    ierr = VecNorm(GhostDofs,NORM_2,&nrm2_dofs); CHKERRQ(ierr);
    ierr = VecNorm(Ghost_b,NORM_2,&nrm2_b);  CHKERRQ(ierr);
    ierr = PetscPrintf(pcomm->comm(),"driver_setup_arc_length nrm2_dofs=%6.4e nrm2_b=%6.4e\n",nrm2_dofs,nrm2_b); CHKERRQ(ierr);
    //printf("dofs %6.4e\n",nrm2_dofs);
    //printf("b %6.4e\n",nrm2_b);
    if(nrm2_b==0) {
      ret_val = 1;
    } else {
      double scale_b = 1./nrm2_b;
      ierr = VecScale(Ghost_b,scale_b); CHKERRQ(ierr);
    }
  } else {
    make_bc_arc_length_b_vector(*(ctx->moab),ctx->meshset,Ghost_b);
    VecGhostUpdateBegin(Ghost_b,INSERT_VALUES,SCATTER_FORWARD);
    VecGhostUpdateEnd(Ghost_b,INSERT_VALUES,SCATTER_FORWARD);
    double nrm2_b;
    ierr = VecNorm(Ghost_b,NORM_2,&nrm2_b);  CHKERRQ(ierr);
    ierr = PetscPrintf(pcomm->comm(),"driver_setup_arc_length nrm2_b=%6.4e\n",nrm2_b); CHKERRQ(ierr);
  }   
  VECDESTROY(Ghost_b);
  if(ctx->is_hierarchical) {
    VECDESTROY(Ghost_b_hierarchical);
  }
  VECDESTROY(GhostDofs);
  ierr = VecAssemblyBegin(ctx->b); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(ctx->b); CHKERRQ(ierr);
  ierr = VecCopy(ctx->b,ctx->d_b); CHKERRQ(ierr);
  VecZeroEntries(ctx->delta);
  Vec Ghost_delta;
  if(ctx->is_hierarchical) {
    make_hierarchical_petsc_createGhostVecArray_formLargeVec(*(ctx->moab),ctx->meshset,ctx->delta,&Ghost_delta,nb_lambda_dofs);
    VecGhostUpdateBegin(Ghost_delta,INSERT_VALUES,SCATTER_FORWARD);
    VecGhostUpdateEnd(Ghost_delta,INSERT_VALUES,SCATTER_FORWARD);
    ierr = VECDESTROY(Ghost_delta); CHKERRQ(ierr);
  }
  make_aux_petsc_createGhostVecArray_formLargeVec(*(ctx->moab),ctx->meshset,ctx->delta,&Ghost_delta,nb_lambda_dofs);
  make_bc_arc_length_direction_vector(*(ctx->moab),ctx->meshset,ctx->arc_length_direction,Ghost_delta);
  VecGhostUpdateBegin(Ghost_delta,INSERT_VALUES,SCATTER_FORWARD);
  VecGhostUpdateEnd(Ghost_delta,INSERT_VALUES,SCATTER_FORWARD);
  ierr = VECDESTROY(Ghost_delta); CHKERRQ(ierr);
  ierr = VecScale(ctx->delta,ctx->eps); CHKERRQ(ierr);
  ierr = VecAssemblyBegin(ctx->delta); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(ctx->delta); CHKERRQ(ierr);
  return ret_val;
}
PetscErrorCode driver_solve_snes(SNES snes,snes_ctx *ctx,Vec DfromArray) {
  ParallelComm* pcomm = ParallelComm::get_pcomm(ctx->moab,MYPCOMM_INDEX);
  make_bc_remove_normal_tag_form_skin_faces(*(ctx->moab),0);
  make_bc_mark_skin_if_one_neighbour_elem_faces(*(ctx->moab),0,PETSC_NULL,1);
  make_insert_crack_front_edges(*(ctx->moab),0);
  PetscLogDouble t1,t2;
  PetscLogDouble v1,v2;
  ierr = PetscGetTime(&v1);CHKERRQ(ierr);
  ierr = PetscGetCPUTime(&t1);CHKERRQ(ierr);
  ierr = SNESSolve(snes,PETSC_NULL,DfromArray); CHKERRQ(ierr);
  ierr = PetscGetTime(&v2);CHKERRQ(ierr);
  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
  int its;
  ierr = SNESGetIterationNumber(snes,&its); CHKERRQ(ierr);
  ierr = PetscPrintf(pcomm->comm(),"number of Newton iterations = %D\n",its); CHKERRQ(ierr);
  PetscSynchronizedPrintf(pcomm->comm(),"SNESSolve:: Rank %d Solution Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
  PetscSynchronizedFlush(pcomm->comm());
  Vec D,D_hierarchical = PETSC_NULL;
  int nb_lambda_dofs;
  make_aux_nb_Lambda_dofs(*(ctx->moab),ctx->meshset,&nb_lambda_dofs);
  make_aux_petsc_createGhostVecArray_formLargeVec(*(ctx->moab),ctx->meshset,DfromArray,&D,nb_lambda_dofs); 
  if(ctx->is_hierarchical) {
    make_hierarchical_petsc_createGhostVecArray_formLargeVec(*(ctx->moab),ctx->meshset,DfromArray,&D_hierarchical,nb_lambda_dofs); 
    VecGhostUpdateBegin(D_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
    VecGhostUpdateEnd(D_hierarchical,INSERT_VALUES,SCATTER_FORWARD); 
  } else {
    VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD);
    VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD);  
  }
  double Psi;
  make_hierarchical_calulate_elastic_energy_and_save_elem_data_on_tags_hierarchical(ctx->lambda,ctx->mu,*(ctx->moab),ctx->meshset,D,D_hierarchical,&Psi);
  ierr = PetscPrintf(pcomm->comm(),"Psi = %10.9e\n",Psi); CHKERRQ(ierr);
  VECDESTROY(D);
  if(ctx->is_hierarchical) {
    VECDESTROY(D_hierarchical); 
  }
  return 0;
}
PetscErrorCode driver_snes_destroy(SNES snes,snes_ctx *ctx) {
  ierr = VECDESTROY(ctx->b); CHKERRQ(ierr);
  ierr = VECDESTROY(ctx->d_b); CHKERRQ(ierr);
  ierr = VECDESTROY(ctx->delta); CHKERRQ(ierr);
  ierr = VECDESTROY(ctx->dF_dLambda); CHKERRQ(ierr);
  ierr = VECDESTROY(ctx->dD_dF_dLambda); CHKERRQ(ierr);
  ierr = SNESDESTROY(snes); CHKERRQ(ierr);
  #ifdef PETSC_3_3_
  ierr = SNESLINESEARCHDESTROY(ctx->linesearch); CHKERRQ(ierr);
  #endif
  ierr = PCDESTROY(ctx->pc); CHKERRQ(ierr);
  ierr = PCDESTROY(ctx->pc_uzawa); CHKERRQ(ierr);
  return 0;
}
void driver_data_on_mesh(Interface& moab,EntityHandle meshset,snes_ctx *ctx,Vec D,Vec D_hierarchical,Vec F) {
  if(D != PETSC_NULL) make_aux_save_data_on_mesh_nodes(moab,meshset,D);
  if(D_hierarchical != PETSC_NULL) make_aux_save_data_on_tags_for_hierarchical_dofs(moab,meshset,D_hierarchical);
  if(F != PETSC_NULL) make_aux_save_forces_in_tags(moab,meshset,F);
}
void driver_calulate_forces(SNES snes,Interface& moab,EntityHandle meshset,snes_ctx *ctx,Vec DfromArray,Vec FfromArray) {
  Vec D,D_hierarchical = PETSC_NULL;
  Vec F,F_hierarchical = PETSC_NULL;
  VecZeroEntries(FfromArray);
  VecAssemblyBegin(FfromArray);
  VecAssemblyEnd(FfromArray);
  int nb_lambda_dofs;
  make_aux_nb_Lambda_dofs(moab,meshset,&nb_lambda_dofs);
  make_aux_petsc_createGhostVecArray_formLargeVec(moab,meshset,DfromArray,&D,nb_lambda_dofs);
  make_aux_petsc_createGhostVecArray_formLargeVec(moab,meshset,FfromArray,&F,nb_lambda_dofs);
  if( ctx->is_hierarchical) {
    make_hierarchical_petsc_createGhostVecArray_formLargeVec(moab,meshset,DfromArray,&D_hierarchical,nb_lambda_dofs);
    make_hierarchical_petsc_createGhostVecArray_formLargeVec(moab,meshset,FfromArray,&F_hierarchical,nb_lambda_dofs);
  }
  set_type_of_analysis(spatail_analysis|material_analysis);
  unset_material_m_blocade();
  make_bc_remove_normal_tag_form_skin_faces(moab,0);
  make_bc_displacements_on_all_nodes(moab,0,PETSC_NULL,bc_free_all,0);
  if(ctx->is_hierarchical) {
    VecGhostUpdateBegin(F_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
    VecGhostUpdateEnd(F_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
    make_hierarchical_elastic_fint_MBTET_hierarchical(
    ctx->lambda,ctx->mu,ctx->alpha,ctx->gamma,ctx->ave_Vol,moab,meshset,
    D,NULL,F,NULL,D_hierarchical,NULL,F_hierarchical,NULL,FREE_TO_MOVE);
    VecGhostUpdateBegin(F_hierarchical,ADD_VALUES,SCATTER_REVERSE);
    VecGhostUpdateEnd(F_hierarchical,ADD_VALUES,SCATTER_REVERSE);
    VecGhostUpdateBegin(F_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
    VecGhostUpdateEnd(F_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
  } else {
    VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD);
    VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD);
    make_elastic_fint_MBTET(
    ctx->lambda,ctx->mu,ctx->alpha,ctx->gamma,ctx->ave_Vol,moab,meshset,D,NULL,F,NULL,FREE_TO_MOVE);
    VecGhostUpdateBegin(F,ADD_VALUES,SCATTER_REVERSE);
    VecGhostUpdateEnd(F,ADD_VALUES,SCATTER_REVERSE);
    VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD);
    VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD);
  }
  driver_data_on_mesh(moab,meshset,ctx,PETSC_NULL,PETSC_NULL,F);
  make_aux_parallel_exchange_elems_data(moab,meshset);
  VECDESTROY(D);
  VECDESTROY(F);
  if( ctx->is_hierarchical) {
    VECDESTROY(F_hierarchical);
  }
}
void driver_calulate_griffith_forces(SNES snes,Interface& moab,EntityHandle meshset,snes_ctx *ctx,Vec DfromArray,Vec FfromArray) {
  Vec D;
  int nb_lambda_dofs;
  make_aux_nb_Lambda_dofs(moab,meshset,&nb_lambda_dofs);
  make_aux_petsc_createGhostVecArray_formLargeVec(moab,meshset,DfromArray,&D,nb_lambda_dofs);
  driver_calulate_forces(snes,moab,meshset,ctx,DfromArray,FfromArray);
  make_bc_remove_normal_tag_form_skin_faces(moab,0);
  make_bc_displacements_on_all_nodes(moab,0,PETSC_NULL,bc_free_all,0);
  make_bc_mark_skin_if_one_neighbour_elem_faces(moab,0,PETSC_NULL,1);
  make_area_calculate_clear_data();
  make_area_calculate_values(ctx->eps,*(ctx->moab),ctx->meshset,D,MATERIAL_M,true);
  Vec tmpFfromArray,tmpF,tmpF_hierarchical;
  VecDuplicate(FfromArray,&tmpFfromArray);
  make_aux_petsc_createGhostVecArray_formLargeVec(moab,meshset,tmpFfromArray,&tmpF,nb_lambda_dofs);
  if(ctx->is_hierarchical) {
    make_hierarchical_petsc_createGhostVecArray_formLargeVec(moab,meshset,tmpFfromArray,&tmpF_hierarchical,nb_lambda_dofs);
    mult_area_mult_add_shellMat_CT_AREA_by_scalar(tmpF_hierarchical,ctx->gc);
    VecGhostUpdateBegin(tmpF_hierarchical,ADD_VALUES,SCATTER_REVERSE);
    VecGhostUpdateEnd(tmpF_hierarchical,ADD_VALUES,SCATTER_REVERSE);
    VecGhostUpdateBegin(tmpF_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
    VecGhostUpdateEnd(tmpF_hierarchical,INSERT_VALUES,SCATTER_FORWARD);
  } else {
    mult_area_mult_add_shellMat_CT_AREA_by_scalar(tmpF,ctx->gc);
    VecGhostUpdateBegin(tmpF,ADD_VALUES,SCATTER_REVERSE);
    VecGhostUpdateEnd(tmpF,ADD_VALUES,SCATTER_REVERSE);
    VecGhostUpdateBegin(tmpF,INSERT_VALUES,SCATTER_FORWARD);
    VecGhostUpdateEnd(tmpF,INSERT_VALUES,SCATTER_FORWARD);
  }
  mult_area_mult_add_shellMat_CT_AREA_by_scalar_to_tag(moab,meshset,tmpF);
  VECDESTROY(tmpFfromArray);
  VECDESTROY(tmpF);
  make_insert_crack_calculate_Griffith_force(moab,meshset,ctx->gc,
    &ctx->min_g,&ctx->max_g,&ctx->ave_g,&ctx->stdev_g,
    &ctx->min_dot,&ctx->max_dot,&ctx->ave_dot,&ctx->stdev_dot,true,ctx->gc_threshold);
  VECDESTROY(D);
}
void driver_data_on_mesh_arc_length(Interface& moab,EntityHandle meshset,snes_ctx *ctx,Vec DfromArray) {
  Vec GhostDofs;
  int nb_lambda_dofs;
  make_aux_nb_Lambda_dofs(moab,meshset,&nb_lambda_dofs);
  make_aux_petsc_createGhostVecArray_formLargeVec(moab,meshset,DfromArray,&GhostDofs,nb_lambda_dofs);
  make_aux_save_data_on_mesh_nodes_arc_length(moab,meshset,GhostDofs);
  Vec GhostDofs_hierarchical;
  if(ctx->is_hierarchical) {
    make_hierarchical_petsc_createGhostVecArray_formLargeVec(moab,meshset,DfromArray,&GhostDofs_hierarchical,nb_lambda_dofs);
    make_aux_save_data_on_tags_for_hierarchical_dofs_arc_length(moab,meshset,GhostDofs_hierarchical);
  } 
  VECDESTROY(GhostDofs);
  VECDESTROY(GhostDofs_hierarchical);
}
void driver_set_data_from_mesh_arc_length(Interface& moab,EntityHandle meshset,snes_ctx *ctx,Vec DfromArray) {
  Vec GhostDofs;
  int nb_lambda_dofs;
  make_aux_nb_Lambda_dofs(moab,meshset,&nb_lambda_dofs);
  make_aux_petsc_createGhostVecArray_formLargeVec(moab,meshset,DfromArray,&GhostDofs,nb_lambda_dofs);
  make_aux_set_DofsVec_form_tags_values_arc_length(moab,meshset,GhostDofs,set_space|set_mat);
  Vec GhostDofs_hierarchical;
  if(ctx->is_hierarchical) {
    make_hierarchical_petsc_createGhostVecArray_formLargeVec(moab,meshset,DfromArray,&GhostDofs_hierarchical,nb_lambda_dofs);
    make_aux_set_DofsVec_hierarchical_from_tags_arc_length(moab,meshset,GhostDofs_hierarchical);
  } 
  VECDESTROY(GhostDofs);
  VECDESTROY(GhostDofs_hierarchical);
}
double driver_solution_with_alpha_r_adaptation(SNES snes,snes_ctx *ctx,Mat M,Vec FfromArray,Vec DfromArray) {
  ParallelComm* pcomm = ParallelComm::get_pcomm(ctx->moab,MYPCOMM_INDEX);
  int nb_lambda_dofs;
  make_aux_nb_Lambda_dofs(*(ctx->moab),ctx->meshset,&nb_lambda_dofs);
  Vec D,D_hierarchical = PETSC_NULL,Lambda;
  make_aux_petsc_createGhostVecArray_formLargeVec(*(ctx->moab),ctx->meshset,DfromArray,&D,nb_lambda_dofs);
  make_aux_petsc_createLambdaVecArray_formLargeVec(*(ctx->moab),ctx->meshset,DfromArray,&Lambda,nb_lambda_dofs);
  if(ctx->is_hierarchical) {
    make_hierarchical_petsc_createGhostVecArray_formLargeVec(*(ctx->moab),ctx->meshset,DfromArray,&D_hierarchical,nb_lambda_dofs);
  }
  #ifdef PETSC_3_3_
  ierr = PCReset(ctx->pc); CHKERRQ(ierr);
  #endif
  //ierr = SNESReset(snes); CHKERRQ(ierr);
  ierr = SNESGetKSP(snes,&(ctx->ksp)); CHKERRQ(ierr);
  //ierr = KSPReset(ctx->ksp); CHKERRQ(ierr);
  PC pc;
  ierr = KSPGetPC(ctx->ksp,&pc); CHKERRQ(ierr);
  //ierr = MatSetUnfactored(ctx->B); CHKERRQ(ierr);
  //ierr = MatSetUnfactored(ctx->Bprec); CHKERRQ(ierr);
  MatStructure flag = SAME_NONZERO_PATTERN;
  ierr = PCSetOperators(pc,M,M,flag); CHKERRQ(ierr);
  double alpha0 = ctx->alpha;
  double kappa0 = ctx->kappa;
  double arc_lambda;
  int max_it = 6;
  for(int it = 0;max_it!=0;max_it--,it++) {
    //PetscAttachDebugger();
    PetscPrintf(pcomm->comm(),"R-ADAPTATION it = %d alpha = %6.4e \n",max_it,ctx->alpha); 
    make_aux_set_alpha(*(ctx->moab),0,ctx->alpha,16,8);
    make_aux_get_arc_lambda(*(ctx->moab),ctx->meshset,Lambda,&arc_lambda);
    driver_data_on_mesh(*(ctx->moab),ctx->meshset,ctx,D,D_hierarchical,PETSC_NULL);
    //make_aux_set_mesh_coords_fromTags(*(ctx->moab),0,0);
    //scalinig
    rhs_snes(snes,DfromArray,FfromArray,ctx);
    lhs_snes(snes,DfromArray,&M,&M,&flag,ctx);
    ierr = PCSetUp(pc); CHKERRQ(ierr);
    double dot_u_lambda_F_lambda;
    ierr = VecDot(ctx->dF_dLambda,ctx->dD_dF_dLambda,&dot_u_lambda_F_lambda); CHKERRQ(ierr);
    double dot_u_lambda_b;
    ierr = VecDot(ctx->d_b,ctx->dD_dF_dLambda,&dot_u_lambda_b); CHKERRQ(ierr);
    if(fabs(dot_u_lambda_b)>0) {
      double scale = dot_u_lambda_F_lambda/dot_u_lambda_b;
      PetscPrintf(pcomm->comm(),"    scale dot_u_lambda_F_lambda = %6.4e dot_u_lambda_b = %6.4e scale = %6.4e\n",dot_u_lambda_F_lambda,dot_u_lambda_b,scale); 
      ierr = VecScale(ctx->b,scale); CHKERRQ(ierr);
      ierr = VecScale(ctx->d_b,scale); CHKERRQ(ierr);
      ctx->arc_diagonal *= scale;
      ctx->arc_length_value *= scale;
    }
    //end<-scaling
    driver_solve_snes(snes,ctx,DfromArray);
    SNESConvergedReason reason;
    SNESGetConvergedReason(snes,&reason);
    if(reason<=0) {
      PetscPrintf(pcomm->comm(),"Set dof from tags form it %d\n",it-1); 
      make_aux_set_DofsVec_form_tags_values(*(ctx->moab),ctx->meshset,D,set_space|set_mat);
      if(D_hierarchical!=PETSC_NULL) {
	make_aux_set_DofsVec_hierarchical_from_tags(*(ctx->moab),ctx->meshset,D_hierarchical);
      }
      make_aux_set_arc_lambda(*(ctx->moab),ctx->meshset,Lambda,arc_lambda);
      break;
    }
    if(ctx->alpha<=ctx->alpha_min) break;
    int iter; 
    SNESGetIterationNumber(snes,&iter);
    if(iter<=4) ctx->kappa *= ctx->c_inc;
    else ctx->kappa *= ctx->c_red;
    ctx->alpha *= pow(ctx->beta,ctx->kappa); 
    if(ctx->alpha<ctx->alpha_min) ctx->alpha = ctx->alpha_min;
  }
  double alpha_return = ctx->alpha;
  ctx->alpha = alpha0;
  ctx->kappa = kappa0;
  VECDESTROY(D);
  if(ctx->is_hierarchical) {
    VECDESTROY(D_hierarchical);
  }
  VECDESTROY(Lambda);
  return alpha_return;
}
bc_type bc_free_all(int boundcond,double *coord,double *val,bc_type bc) {
  bc_type my_bc = FREE_TO_MOVE;
  if(bc & CRACKED_NODE) {
    my_bc |= CRACKED_NODE;
  }
  if(bc & MATERIAL_F) {
    my_bc |= MATERIAL_F;
  }
  if(bc & MATERIAL_E) {
    my_bc |= MATERIAL_E | ( bc & (MATERIAL_M|MATERIAL_F) );
  } 
  if(bc & ZIPPED_ENTITY) {
    my_bc |= ZIPPED_ENTITY;
  }
  if(bc & UNZIPPED_ENTITY) {
    my_bc |= UNZIPPED_ENTITY;
  }
  return my_bc;
}
bc_type bc_block_material(int boundcond,double *coord,double *val,bc_type bc) {
  return MATERIAL_X|MATERIAL_Y|MATERIAL_Z;
}
bc_type bc_block_space(int boundcond,double *coord,double *val,bc_type bc) {
  return FIXED_X|FIXED_Y|FIXED_Z;
}
bc_type bc_block_material_not_at_fornt(int boundcond,double *coord,double *val,bc_type bc) {
  if(bc&MATERIAL_E) return bc;
  return bc|MATERIAL_X|MATERIAL_Y|MATERIAL_Z;
}


